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ABSTRACT 



Various computer-oriented numerical schemes are 
available to evaluate the velocity potentials for 
inviscid fluid flow past submerged objects. The use of 
a surface distribution of sources of uniform strength 
over panels representing the object’s surface requires 
a relatively fine grid to obtain accurate results. A 
new variation of this scheme using distributed sources 
of linearly varying strengths over triangular panels 
appears to. reduce the number of panels needed for good 
results, thereby reducing computer time and storage 
requirements. In the present study, this technique is 
applied to the problem of determining fluid forces and 
dynamic response resulting from interaction between an 
arbitrarily shaped object and surface waves. 
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I. 



INTROD UCTI ON 



Ocean waves interacting with large marine structures 
produce substantial hydrodynamic forces. These forces are 
experienced by floating objects such as ships or moored 
platforms as well as submerged bodies near the free surface, 
such- as submarines at periscope depth or bottom mounted 
caissons in shallow water. 

It is desirable from an engineering standpoint to be 
able to determine the wave forces and resulting dynamic 
response for any large ocean structure under design. Until 
recently, however, such determinations were not possible in 
the general case. Application of the well-known Mcrison 
equation to calculate wave forces is limited to cases in 
which structure dimensions are small compared to wave 
length, as in the case of a small diameter pile. 

The usual analytical approach in dealing with the 
hydrodynamic interaction of large bodies with waves is based 
on a mathematical modeling of the fluid flow using potential 
flow theory, which assumes the fluid to be inviscid and 
irrotational . While such an assumption is unrealistic for 
steady flow past bluff bodies (because of boundary layer 
separation effects) , experiments indicate that it yields 
valid results for the case of unsteady flow past an object 
as produced by surface waves provided wave amplitude is 
small compared to the characteristic dimensions of the 
object . 

Numerical schemes based on the use of the Green's 
function (cr source distribution) method which utilize 
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high-speed computer capabilities have had considerable 
success in solving problems of wave/body interactions. 
Garrison, et. al. [Refs. 1-5] have developed numerical 
procedures for calculation of wave interaction with large, 
three-dimecsicnal, fixed and floating bodies. These 
procedures, which represent the immersed surface by a grid 
of distributed source panels of uniform strengths, have 
proven to be powerful in practical application as 
demonstrated in Refs, 6-9. 

The most severe limitation of the above mentioned 
procedures is the rather large amount of computer (C.P.U.) 
time reguired to obtain an accurate solution . Thus, any 
improvements to the method which would improve convergence 
of the solution would be of considerable value. 

A natural extension of the technique based on a uniform 
distribution of sources over each grid panel is to allcw the 
source strength to vary linearly over the extent of each 
panel. Since this is not convenient when using quadrilateral 
panels, it is appropriate to divide the source surface into 
a grid of triangular panels. The linearly varying source 
strength may then be represented in terms of the values of 
the source strength at the corners of the triangles. 

Webster [Ref. 10] has developed such a scheme for the 
case of uniform flow of an unbounded fluid past a 
three-dimensional body. The body was represented by use of 
the three-dimensional source potential of the form 1/R where 
R denotes the distance from the source. In the present study 
this general method is extended to the problem of 
representing a body of arbitrary shape located in fluid of 
finite depth in the presence of a free surface. 
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The tasic approach used in solving the problem of wave 
interaction with large floating objects is to decompose it 
into component parts which are more simply analyzed, and 
then ouperimpcse the results to yield the overall solution. 
A group of boundary-value problems are generated in terms of 
velocity potential; one problem is solved for each of the 
six possible degrees of freedom of body motion, and one for 
the scattering of the incident wave by the body. The 
solutions are then formulated in terms of the Green’s 
function. The resulting integral equation is then converted 
to a system of linear equations which are subsequently 
solved by ccuputer. 
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II, iOIMLATIM 0.1 the wave interaction problem 



A. restrictions and assumptions 



All of the usual simplifying assumptions incorporated in 
linear gravity wave theory are adopted in the present study. 
In summary, these are: 

m The flow is considered to be incom pres sible, 
irrotaticnal , and inviscid. 

■ The wave amplitude to wave length ratio is 

small . 



■ The pressure is uniform above the free surface, 
and the density is uniform throughout the fluid. 

■ The bottom boundary is represented by a 

horizontal, impermeable plane. 

■ The velocity-squared term in Bernoulli's 

equation is neglected. 

■ The waves are regular and periodic. 



In addition, it is assumed that 

■ The object has a smooth, rigid, and impermeable 

surface . 

The wave amplitude is small compared to the 
object's dimensions. 



B. PROBLEM EESCRIPTION 



The general description of the problem is shown in Fig. 
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1, in which two coordinate systems are defined. The inertial 
coordinate system, positioned at the mean free surface, has 
directions defined by x, y and z, while the coordinate 
system positioned on the object has directions indicated by 
X*, y* and z*. The origin of the body coordinates is taken 
to be at the center of gravity of the body, and the body is 
submerged an arbitrary distance, d, below the mean free 
surface. The incident wave is assumed to propagate in the 
direction defined by the angle ^ indicated in Fig. 1. 

The analysis of the floating object located in a train 
of regular waves is decomposed into seven problems which are 
uncoupled and can be solved independently and superimposed 
to obtain the complete solution. 

The first six problems are identical in fori, and 
consist of determining the fluid forces due to oscillation 
of the object in the six degrees of freedom, in an otherwise 
still fluid. The six modes of oscillation are all considered 
to take place with angular frequency, C*=Z7r/T, which matches 
the frequency of the incident waves. These harmonic motions 
of the object may be expressed as 



o o o 

in which X =a 0 when k=4, 5, 6, and where 0 denotes the 
k k k 

angular displacement about the x», y* and z* axes, 

respectively. The subscripts k=1, 2, 3 correspond to surge, 
heave, and sway modes of motion, respectively, while the 
subscripts k=4, 5, 6 correspond to the rotational modes of 

roll, yaw, and pitch, respectively, as depicted in Fig. 2. 
o 

The parameter denotes the complex amplitude of the linear 
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modes of motion, while the term a 9 denotes rhe complex 

k 

amplitude of the rotational degrees of freedom. 



The seventh problem to be considered is the interaction 
of waves with the object held fixed in position. This 
interaction cf an incident wave with the fixed body results 
in "scattered” waves which radiate from the body. Prom the 
solution of this problem the excitation forces and moments 
acting on the bo-dy are determined. 



C. 2QUATICNS OP 210TI0N 



For a flcatinc body with linear elastic constraints the 
equation of motion for the six degrees of freedom is given 
in fief. 1 as 




where , and CT =21X/t denotes the frequency of the 

incident wave and resulting motion. (In Eq. 2 and all 
equations which follow, the use of the lower case letter, i, 
as a subscript denotes an index. Appearance of i elsewhere 
in an equation denotes the imaginary number i= 1 •) 

surface elevation of the incident wave at the origin is 
given by 



cos (crt) 



in which is a real number denoting the half-amplitude of 
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th€ incident wave. All phase angles are referred to the 
wave crest with positive values indicating lag. 

In Eg. 2, m denotes the dimensionless body mass tensor 

ij 

and is defined by 






(4a) 






(4b) 



in which m denotes the mass of the body, p denotes the fluid 
density and a denotes the characteristic length scale of the 
body. Similaily, 



where the mcnents of inertia are defined typically as 





' x' y' Jrfj 


(5a) 




y'z.' Jfn 

1 


(5b) 


in which the integration 


extends over the complete 


mass of 



the body denoted by m. 

The excitation force and moment coefficients are defined 
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as 




F,- (ma>0 e ' 

aP V 



J - 




F- Cmax)e'^J 



J = -i', 6 



(6) 

(7) 



where ij denotes the phase shift of the force with respect to 
the crest cf the incident wave. A positive value of Sj 
denotes a lag. The term F (max), (j=1r 2, ...6), denotes 

j 

the maximum value of the force (j=lr 2, 3) and moment (j=4. 



5, 6) acting cn the body. 



In the case of the forces (and moments) due to the 
motion of the body, it is more common to describe the 
hydrodynamic force in terras of the added mass tensor, ia 

ij 

and damping tensor N . These parameters are defined as 
i 3 



— M « r k I , - 0 * ^ I i 






- - i N;.' = 



where F (max) denotes the amplitude of the force or 



(8) 

(9) 

moment 
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in the i-th direction resulting from the motion of the 
floating tody in the j-th degree of freedom. 



The forces defined in Egs. 6-9 are given in terms 
surface integrals cf the hydrodynamic pressure as 



of 






o — I ^ ^ J / * . Cj 



( 10 ) 
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where F denotes the i-th component of wave excitation force 
i 



or moment and F denotes the i-th 

ij 

moment arising from the j-th component 


component of force or 
of body motion. The 


functions h are defined as 
i 






(12a) 


= fly 


(12b) 


hj = 


(12c) 


= Tdi-y) -z ny 


(12d) 


= z " X 1^2. 


(128) 


= xny -(dl-t-y) 


(12f) 



in which n , n and n denote the three components of the 
X y z 

outward unit normal vector on the surface of the body. 



The parameters P , (j=1/ 2, ...6), denote the pressure 

j 

on the immersed surface associated with motion in the six 

degrees of freedom and P denotes the pressure associated 
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with the incident wave (0) and scattering (7) cf the 

incident wave. Finally, in Eg. 2 the dimensionless 
coefficients K denote the hydrostatic restoration force 

ij 

(or moment) tensor and are defined as 









s = 


(13a) 




Kz 




^ z' tly ds 


(13b) 








'n^-z*^ny)ds 


(13c) 








/'z' riy (As 


(13d) 








I x' Hy (As 


(13e) 








n,- x'^ny) ds 


(13f) 


denotes the waterline area. 





w 

Mooring line reactions depend, in general, on the 
mooring line configuration, weight, shape (catenary) , 
hydrodynamic and elastic properties. However, it is 
generally permissible to disregard dynamic effects and 
approximate the reactions by a linear relationship with the 
six components of body displacements as indicated by Eg. 2. 
The term K* in Eg. 2 represents the dimensionless force 

i j 

(or moment) produced in the (negative) i-direction by a unit 

displacement of the body in the j-direction. Given the 
mooring line configuration and elastic properties, the 
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spring constant tensor can be evaluated. 



D. REPRESENTATION BY VELOCITY POTENTIALS 

The fluid motion in each of the seven flow situations 
previously described can be represented by a velocity 
potential. Specification of the velocity potential of a 
flew completely describes the flow, since fluid velocities 
at any point are obtainable through the relation 



where q(x,y,z,t) represents the velocity vector, ^(x,y,z,t) 
represents the time dependent velocity potential, and 

denotes the gradient of the potential. Thus, for the 
six modes of oscillation of the body in still water, the 
fluid velocity vector is given by 



of the total potential induced by the k-th mode of 
oscillation. The velocity potential associated with the 
seventh problem, waves interacting with the fixed body, is 
actually composed of two components. One component is the 
potential due to the linear incident wave, the space 
dependent part of which is given by 



<j_(x,y,z,i) = V^U,y,z,t) 



( 14 ) 





(cosh k(hi-‘r) \ o ^l^'Xcos^r^k:^.sin^l (16) 
cr V cosh kh I 
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Q 

in which ^ =H/2 denotes wave amplitude (half-amplitude) and 
]c= ZTC/l denotes wave number. 

The second component of the wave interaction potential 
is due to the presence of waves whose characteristics have 
been altered by interacting with the object. This potential 
is referred to as the scattering potential and is denoted by 
the subscript (7) . The total potential for the motion of the 
body in regular waves is expressed through the superposition 

§(x>y,z,t)= Z § 

• 1^0 



E. THE ECUNEARY-VALOE PROBLEM 

Each cf the velocity potentials ( <i , k=1,2...7) must 

k 

satisfy the continuity equation, which in the case of 
irrctaticnal flow takes the form 

V^4>^(x,y,z)r 0 , k=o^l.,.,7 (18) 

In additicn, the potentials must also satisfy the kinematic 
bottom boundary condition 

^ k=0,\...7 (19) 

where h denotes the mean water depth, as well as the 
linearized dynamic and kinematic free surface boundary 
conditions expressed by 
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k - 7 



( 20 ) 



|l(c(X,0,z) _ tf^4>^(x,o,z) = o , 

where g (ienotes the acceleration of gravity. 



While the free surface and bottom boundary conditions 
are exactly the same for all of the potentials, the 
kinematic boundary condition applied on the surface of the 
object is different for each case. This boundary condition 
results from the assumption that the surface of the object 
is rigid and impermeable; it simply imposes the condition 
that there be no normal component of fluid velocity relative 
to the surface at any point on the surface, and is expressed 
mathematically as follows: 



For 



^3 

•Sn 




7>y\ 






oscillation of the body in rhe still fluid, 

= 9» = 

= 02,= -'•O'JC* riy 

-i-o'®‘TC(cl+-y)n2: + Zr\y'\ 

- -tcrej Czrix- 

= = -Lcre^ Cxnj, -(d+-pnxl 



(21a) 
(21b) 
(21c) 
(21d) 
(21e) 
(2 If) 



where 

unit 



n , n and n denote the x, y and z componenxs of the 
X y z 

outward normal vector on the immersed surface which is 
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defined in its mean position by S(x,y,z)=0. 



For wave interaction with the fixed body, 

5^7 _q i.a 5 )®k -i(kxcosJf + kz sia y) 

•Cny $mKCk(V\<-y)]+-i (0;t cosy sin y)* Cosh Ck(h^y)]] 

in which X denotes the angle between the incident wave 
propagaticn direction and the x*-axis. Equations 21 and 22 
are applied only on the immersed surface as defined by 
S (x^y,z) =0. 
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III. ^LOTION OF THE BOU NDA RY- VALOF PROBLEM 



A. GREEN‘3 FUNCTION SOLUTION 



In using Green's function representation of the velocity 
potential, each of the seven unknown potentials is 
considered to be induced by fluid sources which are 
distributed over the surface of the body in a continuous 
manner. The strength of the source at any location on the 
body surface is a function of position and is denoted by 

where indicates a point on the immersed 
surface. (Figure 3 illustrates the concept of the continuous 
source disrribution. ) The velocity potential associated with 
each of the boundary value problems is then given by the 
surface integral 



- iC (23) 

s /r=/,2..,7 

where G(x,y,z; f ^ ) is referred to as the Green's 
funcrion. Its value depends upon the locaxion of the point 
( > ) on the source surface, and the point (x,y,z) in the 

fluid field. The Green's function for the present case, 
which satisfies the bottom and free surface boundary 
conditioDS, is given by 

G- = I/R + j/R' + Cr* (2U) 



in which 
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(25a) 



R' ^ (y 

The G* term is given by Ref. 3 as 

^ l^oO 

^ = 2 RV. 

•^o 

(A' cosh [m (V^h)l co%hLi( (yfh)liJ r. CMr)cJM 

^ S/C7/! (Jlh) - rCOShC^^) 

-h i zn (k^-V^) Cosh[k{^r/j)] cosh [k(v^h)l J^(k c) 
k^h -v"^h + V 



in which 



IT +(Z-^)^‘ <27) 

The term J denotes the Bessel function of the first kind, 
0 

of order zero, and V represents dimensionless wave 

frequency. E.V. indicates the principal value of rhe 
infinite integral. 



An alternate series form of the complete Green *s 
function is given by 



G-= ^7t 

k^h -v^h ^■y 



cosh i;ic(hr7)] OoshLlCr^^k^JLV^r) - i Je(f<^)] 



22 



(28) 



+M 



oo 






cos\M^ (/i-h)Jcc3£//^ r) 



in which Y denotes the Bessel function of the second kind, 
0 

of order zero, K denotes the modified Bessel funcxion of 
0 

the second kind, of order zero, and represents the real 
positive roots of the equation 



tan -rl) = 0 



(29) 



In numerical evaluations either form of the Green's 
function may be used with the exception of the case where r 
approaches zero. Here the series form given in Eg. 28 is 
singular so the integral form given by Eg. 26, which has a 
well-defined 1/R singularity, is used. However, when r is 
not small the series form is generally used for numerical 
calculation because it requires less computer time for 
evaluation. 



The Green's function representation of the potential 
given by Eg. 24 satisfies the laplace equation as well as 
the free surface and bottom boundary conditons for all 
cases. The source strength function unknown 

and is determined by application of the kinematic boundary 



conditions on the body surface (Egs. 21 and 22) . This 
results in the following integral equation: 



^ ^5 = 5^(X,y,z) (30) 

tr=t,z-.,7 

in which denotes the deriviative of the Green's 

function in the outward normal direction, which is given by 
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(30a) 



where n represents the outward unit normal vector at {x,j,z) 
and where {x,Y,z) now represents a control point on 
S(x,y,z) =0. 

For the seven problems, all quantities in Eq. 30 are 

known except for the source strength functions f i ^ , 

k ^ 

(k=1,2, ...7). However, the integral equation cannot be 

solved analytically for an arbitrarily shaped body, and a 

numerical approach must be adopted. The numerical approach 

developed herein involves discretizing both the source 

strength function and the body surface boundary condition 

functions, g , (k= 1 ,2, . . . 7) , 
k 

B. DISCRETIZATION OF THE PROBLEM 



• On if crm Source S tren gth Elem ent Method 

The technique previously used [Ref, 1 ] to discretize 
tne integral equations given by Eq. 30 is to divide the 
source distribution into a grid of N quadrilateral panels, 
and to assume that the source strength has a distinct, 
constant value over each panel. Hence, for any of the seven 
flew problems of interest here, f denotes the strength of 

the source everywhere on panel j, (j=1,2...N), for the k-th 
problem (k=1,2 ...7). The area of panel j is denoted by 
ASj and the induced potential anywhere in the fluid thus 
results frem a summation of the effects of all the source 
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panels. Figure 4 depicts the discretized source surface and 
resulting potential. 



Additionally, the normal boundary condition 

functions g (x,y,z), (k=1,2...7), of Eg. 30 are discretized 

k 

by applying the kinematic body surface boundary condition at 



only M selected points on the surface. In practice, the 
number of boundary condition ("control”) points are chosen 
to be equal (M=N) . For convenience, the control points are 
chosen to lie at the centroids of the source panels. Thus, 
for each discretized integral equation, a system of N linear 
equations in N unknowns results: 




(31) 



in which (x ,y ,z ) represents a control point i on 
i i i 

surface, at which the boundary condition is being 
and g denotes the value of g (x,y,z) at point i. 



31 is expressed in matrix form as 



the body 

ap plied. 
Equation 



= [9k} 



k= \ 



(32) 



where [oc] is an NxN matrix composed of elements 



o< 



u = 5 :| ^ <=^1 

The vectors represent the N unknown source 

strengths and the kinematic boundary condition applied at 
the N ccntrcl points, respectively, for the k-th flow 
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problem. It should be noted that the matrix [ o< ] depends 
only on the body configuration, water depth, and period of 
the motion and is the same for all seven cases. 

Ihe discretized potential which is induced by the N 
source panels is specified by a matrix equation similar to 
Eg. 20: 






The martix [j3] is composed of elements 






(35) 



j-y 'ASj 



As in the case of matrix [cx], matrix [^] is the same fcr all 
seven flew cases. 



Integration of the Green’s function and its 
derivatives in Eqs. 33 and 35 poses no great difficulty. It 
is most convenient to separately integrate the singular 
terms. Special ’ considerations in performing these 
integrations and necessary formulas are provided in Ref. 2. 

A computer solution of Eg. 32 is accomplished using 
matrix inversion and multiplication, yielding the previously 
unknown values of the source strengths. Values of velocity 
potential at the control poinrs (panel centroids) are then 
obtained directly using Eqs. 34 and 35. Fluid velocities at 
the control points are similarly obtained after first 
differentiating the elements of the [|G ] matrix with respect 
to the x,y and z directions. Dynamic pressures are 
calculated from the linearized form of Bernoulli’s equation. 
Finally, the forces acting on the object are determined by 
pressure integration (or more precisely, summation) over the 
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surface area of the object as indicated by Egs. 10 and 11. 



2. M^hod of Linearly Var yi ng S ourc e S trength Ele ment s 

An alternative method of discretizing tne Green’s 
function integral equations, proposed by Webster [ Hef . 10] 
for unifcrni unbounded flow past a submerged body, utilizes 
triangular source panels over which the source strength 
function is considered to vary linearly. In this way, the 
discontinuities or ’’jumps" in source strengths at adjoining 
panel edges (which existed in the previous scheme) are 
eliminated. Also, since source strengths are not held 
constant over the extent of each panel, the discretization 
results in a less crude approximation of the continuous 
strength function. An additional "smoothing" effect on the 
velocities induced near the body is achieved by submerging 
the source panel surface inward a selected distance from the 
actual body surface, although the kinematic body boundary 
condition is still applied at points on the true body 
surface. It is noted that while some smoothing effect can 
be accomplished ty submerging the sources, this is not a 
necessity as long as the singularity is accounted for 
properly . 

Figure 5 depicts the linearly varying source 
strength function for a single triangular panel. The local 
coordinate sjstem for each panel is defined such that two 
triangle corners (denoted by a and d) lie on the local 
x-axis, and the third corner (b) lies on the positive 
y-axis. Appendix A gives the transformation between the 
body coordinate system and the local coordinate system for a 
given parel. 

In panel coordinates, the linearly varying source 
strength function may be written 
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f = -=<+(3? 



(36) 



in which ^ y and 'ZT denote constant coefficients and the 
tildes indicate local coordinate values as defined in Fig. 
5. The values of the source strength at the three corners 
are 



constant coei 



corners: 



P , o) 


(37a) 


P(o,b) 


(37b) 


Pca.o) 


(37c) 


^ into Eg. 36 allows the 
expressed in terms cf the 
i of the strength at the 


-ci) 


(38a) 




(38b) 




(38c) 



Thus, the ccnrribution to the potential at (x,y,z) (3ue to 
the single panel of area is denoted by A^CX/YrZ) and is 
given by: 
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^■^AS 



( 39 ) 



where ^ and denote local panel coordinates and (x^y,z) 
denotes a general point expressed in terms of local 
coordinates. 



Substituting Egs. 38 into Eg. 39 and collecting 
terms allows the contribution to the potential to be 
expressed in terms of the values of the source strength at 
the corners of the panel: 



Ac}>Cx,?,z) («0) 



in which 



A (^.7 ; 1 2) c/s (4 1 b) 



for convenience referred to as 
simply the "corner integrals" associated with corners a,b 
and d, respectively, for a given triangular panel. 



The velocity induced by a single source panel is 
obtained by taking the gradient of Eg. 40: 



A^ (x,y,z) 
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The gradients of the corner integrals are given by 




(43a) 



(43b) 



(43c) 



The evaluation of the three corner integrals (Eg. 



41) and their gradients (Eg. 43) must be accomplished in 
local coordinates for each panel. In performing the 
integrations^ it is most efficient to separate the Green's 
function into its singular and non-singular components, and 
treat the integrals separately, as follows: 



-{( Wj ('/R^dls f ([ W4 ('/R')d s )- f( (G-*) ds 
VfAcii^Vf KV(t/R')is^((vyc^^(&*)ds 

''a 5 

V’fA'tt,) ^(!/R')ds-hjJ v(C>^ds 

f- J' V(>A)Jst 
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in which S , 
a 

are given by 



K t and W are linear functions in 
b d 



f and and 





H7T b(a-J) 


(46a) 




J1 


(46b) 


o 


nn b 






ba— Vj ? - a ^ 
HKh(a-d) 


(46c) 



Appendix B provides formulas developed by Webster 
[Sef. 10] for the closed form evaluation of those integrals 
in Egs. 44 and 45 which involve derivatives of 1/R and 1/R* 
as well as similar formulas developed by Yeung and Bai [Ref. 
11] for integrals of 1/R and 1/R’. The integrands are 
singular, and it is therefore necessary to evaluate the 
integrals analytically. 

Integration of the G* components of Egs. 44 and 45 
is made possible by assuming that the G* function and its 
gradient vary linearly over the triangular panel, (as was 
done for the source strength function) . This is a 
reasonable approximation since is a well-behaved function 
that varies slowly with position on the body surface. The 
wave length of its oscillation is approximately the same as 
the wave length of the incident wave, and this tends to be 
large relative to the panel dimensions in typical 
applicat icns . 

With G* and VG* expressed as linear functions of 
their values at the panel corners, the terms in Egs. 44 and 
45 involving these two functions may be expressed as 
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^ (a.M, G-/, 


(a7a) 




(47b) 




(47c) 




(48a) 


V'/, VCr^'cls ^ rj 


(48 b) 




(48c) 



where I , I , I , 1*^ I*, I’/ represent the integrals cf the 
a b d a b d 

product cf the two linear functions over the triangular 
area. The integrals in Egs, 47 and 48 were evaluated using 
the formulas given in appendix C for integrating the product 
of two linear functions over a triangular area and the 
resulting integrals are given in Appendix D 

The total potential and velocity at any point in the 
fluid is a result of the combined effects of all the source 
panels. For a grid composed of a panels, the potential is 
given by 







and the velocity vector is given by 

k k k 
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in which f rf # f / denote the strengths at the three 
^Ic '^K 

corners of the k-th panel. 

It should be emphasized that the corner integrals 
and their gradients are evaluated in their own individual 
panel coordinate systems. Transformation back to the global 
coordinates, described in Appendix A, must be accomplished 
prior to performing the summations of Eqs. 49 and 50. 

The summations given in Egs. 49 and 50 can be 
re-arranged by reversing the order of the summation. That 
is, if coefficients of the strength at a given node are 
first summed, then the summation may be written 

4>(x,y,z) = fj A(J)j (x,y,z) 

f fj 9'(A4>,(x,y,Z)) 

in which f denotes the source strength at corner node j 

j 

(j = l,2...N) and represents the sum of the corner 

integrals associated with node j, for all panels which share 

node j as a mutual corner. A graphic representation of the 
meaning of A(J). is provided in Pig. 6 for clarity, since it 

is difficult to express mathematically. As indicated in 
Fig. 6 A(j>.(x#y,z) denotes the contribution to the 

potential at (x,y,z) from the panels surrounding node j. 

Each of the surrounding panels may be considered to have 
unit strength at the corner coincident with node j and zero 
at the other two corners as depicted in Fig. 6. 



(51) 

(52) 
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The kinematic boundary condition at the body surface 
may now he applied at N selected points. It is convenient 
to use points on the surface which lie directly outward from 
the corner nodes of the submerged source surface. In 
practice, it is best to initially define the nodes on the 
actual body surface, and then establish the points on the 
submerged grid at some pre-selected submergence distance 
along the corresponding inward normal vectors. 

Application of the kinematic boundary conditions for 
each of the seven problems of interest results in the 
discretized version of the integral equation; 



r 

^kj an 




, k=i,2..„Tj i = (53) 



where the repeated index j indicates summation. The term 
denotes the sum of the corner integrals for node j, 
evaluated at control point i. The derivative with respect 
to the outward normal vector is obtained by 



(54) 

'b r\ 

where n is the unit outward normal vector at the control 
i 

point i, cn the body surface. 



Equation 53 may be expressed more simply as 



where the NxN matrix, [o< ], is composed of elements 
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( 56 ) 



_ _L_ •& 



Equation 55 is solved by computer, (for ]c=1,2. ..7) 

using matrix inversion and multiplication, yielding the 

unknown source strengths. Substitution of the source 

strengths into Egs. 51 and 52 with evaluation at (x ,y ,z ) 

i i i 

(i = 1 , 2 , . . . N) , gives the values of the potential and the 
fluid velocities at the control points on the body surface. 



Ihe dynamic fluid pressure is determined using the 
linearized form of Bernoulli's equation. For oscillation of 
the body in its six degrees of freedom, the pressures at the 
surface nodes (denoted by i) are given by: 



P|^. = Re- fi/’cr , k=i,2...7 (S7) 

pQ^. = Re Lifer (45oi 1 



where denotes fluid density, and and 

dynamic pressure and potential, respectively, 

(denoted by the subscript i) due to motion in the 
of oscillaticn. 




at 



k-th 



denote 
node i 
mode 



The pressure integrals given by Egs. 10 and 11 
cannot be utilized directly to obtain total forces and 
moments, because only discrete values of pressures are 
available. However, it is consistent with the numerical 
scheme to assume that the pressure varies linearly over each 
triangular panel. Based on this assumption, the pressure at 
any point on a surface panel is expressed as a linear 
function of the panel geometry and the values at the corner 
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nodes, which are given by Egs. 57 and 58. 



The forces and moments are then calculated using 






- f I a hi js 

i-i ''•’as,- 






Fl (t) = ^ ff Pay h; c/s , C: 



I 



(59) 



(60) 






where F denotes the i-th component of wave excitation force 
i 

or moment and F denotes the i-th component of force or 
ik 

moment arising from the k-th component of body motion. The 



summations of Egs. 59 and 60 are carried out over the M 

panels on the actual body surface. The functions h are 

i 

given by Eg. 12 and are functions of position on the 

immersed surface. It appears to be consistent with the 
numerical scheme to approximate these functions as linear 
functions so that they may be defined at interior points by 
their values at the corners of the triangles. Thus the 
integrands of Egs. 59 and 60 contain the products of two 
linear functions. The integrals are evaluated in a 
straight-forward manner using the integration formulas given 
in Appendix C . 



The resulting forces F (i ,k= 1 , 2. . .6 ) and F 

ik i 

(i=1,2...6) are substituted into Egs. 6, 1, 8 and 9 to 
evaluate the excitation force coefficients C ( j= 1 , 2 , . . . 6) , 

j 

and the added mass and damping coefficients, M , and N , 

ij ij 

(i # j= 2 , . . . 6) . Finally, the eguations of motion for the six 
degrees cf freedom (Eg. 2) are applied to determine the 
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ccniplete dynaisic 



response of 



the body. 
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IV. 



NUMERICAL 



A computer code based on the numerical method outlined 
was developed to calculate rhe hydrodynamic coefficients as 
well as the dynamic response of a floating body in waves. 
Using this program, example calculations were made for two 
simple geometric shapes; a vertical circular cylinder which 
extends from the bottom and passes through the free surface, 
and a semi-immersed sphere. Computer produced drawings of 
the grid configurations used are shown in Figs. 7-10. 

The vertical circular cylinder represents the only 
three-dimensional geometry for which a closed-form solution 
exists (of the type of interest in the present study), and 
therefore, it is of interest for making a comparison of 
numerical results computed by the present method. MacCamy 
and Fuchs [Ref. 12] have developed a closed form solution 
for the horizontal force acting on the cylinder and this 
result is plotted in in Fig. 11 in the form of the 
dimensionless amplitude of the force versus the wave length 
parameter, 27Ta/L. 

Corresponding results computed by the triangular panel 
method using grids of two different finenesses are shown on 
Fig. 11 for comparison with the results of MacCamy and 
Fuchs. These results indicate rather rapid convergence of 
the solution, particularily at small values of 2lTa/L. The 
slower convergence (or greater error) at the larger values 
of 2 n a/L apparently indicates the inaccuracies associated 
with the assumption of linearly varying source strength. As 
2 TT a/L becomes large (the wave length becomes small) the 
variation of the part of the Green's function with 
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distance beccmes more rapid. Thus, as 2TTa/L increases, a 
point is reached where the assumption of the linear 
variation of the Green’s function becomes appreciably 
invalid. Of course, as the fineness of the grid is 
increased the numerical results tend to agree with the 
clcsed-fcrm exact results up to higher values of 2Ha/L. 

Havelock [Ref. 13] has presented results for the added 

mass and damping coefficient in heave for a semi-immersed 

sphere and these results are presented in Fig. 12 for 

comparison with numerical results computed by the present 

method. The figure shows the heave added mass and damping 

coefficients, M and N , respectively, for the 

22 22 

semi-immersed sphere along with the numerical results 

computed using grids of two different finenesses. These 
results indicate that the numerical results are tending to 
converge although even for the finer grid the results are 
not yet converged. 

Moreover, given the damping coefficient, the wave force 
coefficient can be computed by use of the Haskind’s 
relations which, in general, relate the force coefficients 
to the £ar-fieid solutions of the corresponding radiation 
problem. For the case of the heaving motion of an 
axi-symmetric body a closed form relationship exists between 
the heave damping coefficient and heave excitation force. 
Thus, Havelock’s results for the damping coefficient in 
heave were used to compute the heave excitation force and 
this result is presented in Fig. 13 for comparison with the 
numerical results. 

The numerical results presented for comparison with 
Havelock’s results in Fig. 13 show a trend toward 
convergence and, in fact, the results corresonding to rhe 
finer grid indicate adequate agreement. 
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Finally, the dynamic response in heave and surge for the 
case of the semi-immersed sphere are shown in Fig. 14. 
These results, showing the dimensionless response in terms 
of the response amplitude to wave amplitude ratio, indicate 
very little difference between the the two grid finenesses. 
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V. 



CONCLUSIONS 



A new numerical procedure based on the use of triangular 
source panels has been developed for computing the 
interaction of fixed or floating bodies with waves. 
Numerical results based on this procedure compare well with 
existing results, and, in general, the procedure appears to 
converge to the correct solution. 

The value of the triangular panel procedure as compared 
to the procedure based on quadrilateral panels of uniform 
strength remains uncertain. This will be determined only 
through experience in application of the two methods. 
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APPENDIX A 



COORDINATE TRANSFORMATION 



In the following, the coordinate transformation is 
developed which relates the local coordinates attached xo a 
particular panel, as indicated in figure la, to the global 
coordinates with origin located at the mean water level. The 
global system is denoted by 0(x,y,z) and the local 
coordinate system is denoted by 0(x,y,z). 




The location of the origin of the local coordinate 
system expressed in global coordinates is given by 
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in which (x ,7. ) , {x /y ,z ) and (x ^y ,z ) are the global 
aaa bbb ddd 

system coordinates of the triangle corners a^b and d, 

respectively^ and where 

o< = - ^a) (x X, A r(Y^- tl) *(l^ -zj)(Zb~2Jl) 

The unit vectors in the x^y and z directions of the 
local coordinate system are given in terms of their 
respective global system components as 

77 ___ 

J = Uy^ I Uy^ J -h 

in which 1, j and IT are the unit vectors in the y and z 
directions^ and 1, j and ic are the unit vectors in the x, y 
and 2 directions. 

The terms u , u ^ ,..u are defined as fellows: 

^X, = (x^-XaVla;! 

(y^-)^d)/i aai 

Uy, - (xt-Xo)/l oWl 



43 



% = (2b-2o)/lobl 

'■'x.j ^¥,“ 

= <Jx, Uy^- 

where 

I +■ (/a - y^) S (Za -Z^f' 

and 



l°^l =V( 5 (=-^o" Kio-yff ^z,-z^f' 

The transformation matrix is then defined as 




(A y ^ X ^ X 
A j ^ Z ^3 

Uy tv ^ 
y y / j. yj 



Transformation between global coordinates and local 
coordinates is accomplished in the following manner: 



r ^ 

X x-x, 

'7 =[T] jy-Xo. 

7) [z-zj 



'x' 

Y 


= [rf < 


X 

y 


+ 1 


'x; 

Vo 


z 

J 




z 

J 




- / 
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r 

in which [T] denotes the transpose of the matrix [T]. 



Similarly/ a vector quantity, F, evaluated in the local 
coordinate system is transformed bacK to global coordinates 
using 




in which F , F and F are the x, y and z components in 
X y z 

local coordinates and F , F and F are the corresponding 

X y z 

components in the global system. 
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APPENDIX B 



EOSJIULAS FOH EVALUATING INTEGRALS INVOLVING 1/R AND 1/R* 



The terms or Eg. 44 which involve 1/R and 1/B* are 
evaluated as follows, using the panel geometry defined by 
Fig. 1b. 
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ff _ / r ^ A 

])as ‘f^bla-cOR ~ zb(a-d)i^-n T- aR^A 

+(b^ ad)Fab/e" - ^ C 2 >>^'- (y^^ad)(z^^^^)]il/e 

+^( 2 ^--^) lab Ciyp-a(y %z^)] X^a ] 

- 7bXir[2-y2T + R^^afF /c^-d/^^/e^ 
+ C 2 .yJ 5 -;a rzV-^)]r^/c ~ ( Y z"^) 
-Lay^'-c/(zV^")]r;b/e 

in which the following definitions are applicable 



r"^ ~^|{.x-a■')^ f y- r £^' 
v'b Vx^ ^(y--b)'^ t Z^ 



c ~-\J i- b^' 



/> — b/( -ray -ah J^'^b^roYy -c/A 

aX -dz-hf 

4 = ^ 4' --dx - by r 6^ 

+ tan-' 



Zb 



= !°S r rx-^)]/£r^ 



J/, ^ ^ rc/a-"4)/rcr, ±^,)1 

- /.^rfera Z4Vre »b t^')] 
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Rod ~ ^ ^ 

'-4 '^t 



i- _ j, ^ 

^dio ~ ^ah 
^ad " d-^ti 



Formulas for evaluating the terms of Eg. 45 which 
involve tie gradients of the 1/R and 1/R* integrals are 
given by Webster [Ref. 10] as 



aS 



^('/r)a = Jyiitii. ^ (r-r) 
A-TTb(x-d) ' I ^bl 



A5 



+ r A—r f T-^ - 1 

L <: 0 -:l ) c.^ J ^oia e(a-d) ^ 

ba-b?"^7 a^n/p) J 5 - ^-ylid^lr _ 
h7iLU-c\) I e-''^ 

+ [ iir-j) - - ^jj ] 

^ * ri;;- '-.<') - 1 



j? ^ hc(a-c'l) di^ bt(p-cl) 
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-i^)- L #' - ^Jt. - %0) J 

i h , . f .(-.-0 









A 

■_/>T_ 






•/" CZ -2^^ /- ^ z -2^^ J 



V-z-r„;] 
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APPENDIX C 



fCBMULAS EOR INTEGRATION OF LINEAR FUNCTIONS OVER TRIANGULAR 

AREAS 



All formulaiS 
positioned as shown 



given are 
in Fig. 1c, 



applicable for a 
with a>d and b>0. 



rr iangle 




Figure 1c- Triangle Geometry 



If f (x,y) and 



as 



f(x,y) = tx + pxi-X/ 

o<' f/3'x 

and the function values ar the triangle corners are denoted 
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by f # f / f / g / g^# g,/ 

a b d a b d 


then: 


^ ^ afd-df^ 

a -d 




I I 

II 


(3'= Sa '‘34 


-(a-d) fij-afj f-d ^ 




b Cd -d) 


bC^-s/) 



and the integrals of f(x,y)ds and f (x,y) g (x,y) ds are given 
by 




3 ' 



and 




— o<cX^ V/| ^(pc<‘ tcX/ 3 ') W2. 



+ (c<'X +oC ^^')W3+(V/3'+ /3 3") +/3/9 'iV5- 



where 

w,- (( 

z 

vi.-{{ xds - 

<. »s ^ 
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ZH 

W5-=ft = (^^v^wW)^ 

/z 

" f{ V 
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APPENDIX D 



INTSGHAIICN FORI^ULAS FOR INTEGRALS INVOLVING AND VG=^ 



The assumption of a linear variation of G* { ; x, y , z) 
over a given triangular panel, and the subsequent 
application of the formula for the integral of the product 
of two linear functions, given in Appendix C, leads to the 
following expressions for Egs, 47 and 48 which involve the 
G* term. (refer to Fig. lb for applicable triangle 
gecmetry . ) 




in which 
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^ . b(xa-d) 




^3 - y , 


B, -Ca-t-d'^b 
ZM 


B.= ii 


B3 - '=’4 


la. = b fad 


hi 


D3 - ^4 




ZH 





The integrations involving the gradients of similarly 
result in the following expressions: 









11^5 V/rfc 



where 



VC, - i. ^ J 2^(d,o') 4.^, ^(a,o) 

' 2 X -ay 3 Z 

— i _ y _r a<s-y^d.o') 

■ 3 X 

VC^= [(a-d')(r 
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a{ ‘ ^ J -^y,o) ^ I, l£%a)] 

' 3X 37 3Z / 



-I- d ( ‘ t j /b 

^X Dy 3Z ^J ^ 

yc^ - [_a.(l^^(d>^) i- Tlf’V.o^ r^^^ljo)] 

^ 3 X J 

-^(l ^ ^ 

■3/ "dz ''J 
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Pig,-1- Problem Definition 
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Fig. 2- Definition of Body Motion 
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y ,n 




z, c 



Fig. 




X,C 



4>=/g(^,n , c)G(^,n,(;;x,y,z)ds 



)- Potential Due to Distributed 
Sources 
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Pig. 4- Uniform Source Strength Discretization 
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f 




X , 



37 ,n 

Fig. 5- Panel of Linearly Varying Source 
Strength 
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• i 




Pig. 



6- Contribution to Potential at Node i 
Prom Panels Surrounding Node j 
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Figure 7 - QUARTER OF A 13 NODE HEMISPHERE 
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Figure 8 - QUARTER OF A 25 NODE HEMISPriSRE 



63 



r 




Figure 9 - QUARTER OF A 48 NODE CIRCULAR CYLINDER 
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Figure 10 - QUABTER OF AN 80 NODE CIRCULAR CYLINDER 
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Fig. 11- Horizontal Force on a Circular Cylinder 




Fig. 12- Heave Added Mass and Damping Coefficients 
for a Floating Hemisphere 
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Pig. 13- Heave Excitation Force for a Floating 
Hemisphere 
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Fig. 14“ Heave and Surge Response for a Semi-immersed Sphere 
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